function pairedAverage(varargin)
% pairedAverage
%
% Take the average of multiple traces. Intended for paired recordings.
% This program assumes that the data structures of all files are the same.
% Notes: data made available as 'currentTrace' variable in the workspace

% TODO:
%   

% editing:
% 
% taro Kiritani aug 2009 --- adapted from qvTrace
% tk jan 2010 -- add arguments. PAIREDAVERAGE(FILENAME,PATHNAME) takes the
% average of the files. FILENAME should be the cell of strings. PATHNAME is
% a string.
% tk jul 2010 -- compatible with multiple recordings (> 2).
% ----------------------------------------------------

defaultDirectory = 'C:\DATA\Taro\CELLS';

data = [];
        
try
    cd(defaultDirectory);
catch
try
    cd('C:\Data');
catch
    disp('qvTrace: unable to find default dir.')
end
end

if nargin == 0
    [filename, pathname] = uigetfile('*.xsg', 'Select an XSG trace','MultiSelect', 'on');
else
    filename = varargin{1};
    pathname = varargin{2};
end


if iscell(filename)
    fileLength = length(filename);
else
    fileLength = 1;
end

if isequal(filename,0) || isequal(pathname,0)
    return
end

for i = 1:fileLength
    fullfilename{i} = fullfile(pathname, filename{i});
    [pathstr{i}, tracename, ext, versn] = fileparts(fullfilename{i});
    data.completename{i} = fullfile(pathstr{i}, [tracename, ext, versn]);
    % open the file, extract the data
    currentTrace{i} = load(data.completename{i}, '-mat');
    currentTrace{i}.userFcn.filename = data.completename{i}; 
end

data.hfig = figure;
set(data.hfig, 'Color', 'w', 'Name', 'pairedAverage', ...
            'NumberTitle', 'off', 'Position', [360    55   560   657]);
        
set(data.hfig, 'UserData', data);



% put the data in the workspace
assignin('base', 'currentTrace', currentTrace);

% % TODO: fix this temporary KLUGE
% if ~isfield(currentTrace{1}.header.acquirer.acquirer, 'acqOnArray')
%     currentTrace{1}.header.acquirer.acquirer.acqOnArray = 0;
% end

% determine what was acq'd and stim'd, and how
NacqEphys = currentTrace{1}.header.ephys.ephys.startButton * ...
    sum(currentTrace{1}.header.ephys.ephys.acqOnArray);
Nother = currentTrace{1}.header.ephys.ephys.startButton * ...
    sum(currentTrace{1}.header.ephys.ephys.stimOnArray);
%     currentTrace{1}.header.acquirer.acquirer.startButton * ...
%     sum(currentTrace{1}.header.acquirer.acquirer.acqOnArray) + ...
%     currentTrace{1}.header.stimulator.stimulator.startButton * ...
%     sum(currentTrace{1}.header.stimulator.stimulator.stimOnArray);

% decide subplot layout
Nsub = 2 * NacqEphys + Nother; % ephys traces count double
plotnum = 0;

% were data obtained by pulseJacker? (If so, different method for determining pulses)
pulseJackerFlag = currentTrace{1}.header.pulseJacker.pulseJacker.enable; 
if pulseJackerFlag; 
    disp(' '); 
    disp('pulseJacker trace'); 
end

% who triggered whom:
disp(' ');
disp(['ephys startButton: ' num2str(currentTrace{1}.header.ephys.ephys.startButton)]);
% disp(['acq startButton: ' num2str(currentTrace{1}.header.acquirer.acquirer.startButton)]);
% disp(['stim startButton: ' num2str(currentTrace{1}.header.stimulator.stimulator.startButton)]);
disp(['ephys externalTrigger: ' num2str(currentTrace{1}.header.ephys.ephys.externalTrigger)]);
% disp(['acq externalTrigger: ' num2str(currentTrace{1}.header.acquirer.acquirer.externalTrigger)]);
% disp(['stim externalTrigger: ' num2str(currentTrace{1}.header.stimulator.stimulator.externalTrigger)]);

% % loop??: yes, if selfTrigger=0 for all started programs (startButton=1)
% loopFlag = 0;
% loopDaqFlag = 0;
% loopCpuFlag = 0;
% loopIterations = currentTrace.header.loopGui.loopGui.iterations;
% loopIterationCounter = currentTrace.header.loopGui.loopGui.iterationCounter;
% numberOfSelfTriggeringPrograms = currentTrace.header.ephys.ephys.startButton * currentTrace.header.ephys.ephys.selfTrigger + ...
% currentTrace.header.acquirer.acquirer.startButton * currentTrace.header.acquirer.acquirer.selfTrigger + ...
% currentTrace.header.stimulator.stimulator.startButton * currentTrace.header.stimulator.stimulator.selfTrigger;
% if numberOfSelfTriggeringPrograms == 0
%     loopFlag = 1;
%     if currentTrace.header.loopGui.loopGui.cpuTiming
%         loopCpuFlag = 1; 
%         loopDaqFlag = 0;
%     else
%         loopCpuFlag = 0; 
%         loopDaqFlag = 1;
%     end
%     loopIterations = currentTrace.header.loopGui.loopGui.iterations;
% end
% disp(['loopFlag: ' num2str(loopFlag) '; loopCpuFlag: ' num2str(loopCpuFlag) '; loopDaqFlag: ' num2str(loopDaqFlag) ...
%     '; iterations: ' num2str(loopIterations) '; counter: ' num2str(loopIterationCounter)]);

% currentTrace.userFcn.fIAnalysis.familyType = '';

% ============ EPHYS TRACES ==================
for channelNum = 1:NacqEphys
    if currentTrace{1}.header.ephys.ephys.startButton && currentTrace{1}.header.ephys.ephys.acqOnArray(1) == 1;

        ephysdata = zeros(fileLength,length(currentTrace{1}.data.ephys.trace_1));
        for i = 1:fileLength
            ephysdata(i,:) = getfield(currentTrace{i}.data.ephys,['trace_',num2str(channelNum)]);
        end

        timedata = 0 : size(ephysdata, 2) - 1;
        timedata = timedata / currentTrace{1}.header.ephys.ephys.sampleRate;
        plotnum = plotnum + 1;
        subplot(Nsub, 1, [plotnum plotnum+1]);
        plotnum = plotnum + 1; % double-high subplot for ephys traces

 
        if size(ephysdata,2) == 10000
           traceBegin = 1;
           traceEnd = 10000;
        end
        if(~exist('traceBegin'))
            traceBegin = inputdlg('beginning point');     
            traceEnd = inputdlg('end point');
            traceBegin = str2num(cell2mat(traceBegin));
            traceEnd = str2num(cell2mat(traceEnd));
        end
        mode = input('Take average (1) or draw each trace (2)? Say 1 or 2: ');

        
        switch mode
            case 1
                meantrace = mean(ephysdata,1);
                plot(timedata(traceBegin:traceEnd), meantrace(traceBegin:traceEnd))
            case 2
                plot(timedata(traceBegin:traceEnd),ephysdata(traceBegin:traceEnd)) % this might be a bug.
        end
        xlabel('Time (seconds)');
        ylabel('Amplitude');
        traceType = '';

    %     if loopFlag
    %         if loopCpuFlag
    %             traceType = [' -- cpu loop, iteration ' ...
    %                 num2str(loopIterationCounter) ' of ' num2str(loopIterations)];
    %             currentTrace.userFcn.fIAnalysis.familyType = 'cpu loop';
    %         elseif loopDaqFlag
    %             traceType = [' -- daq loop, all ' num2str(loopIterations) ' iterations'];
    %             currentTrace.userFcn.fIAnalysis.familyType = 'daq loop';
    %         end
    %     else
    %         currentTrace.userFcn.fIAnalysis.familyType = 'not a loop';
    %     end
    %     htitle = title([pathstr '\' tracename traceType]);
    %     set(htitle, 'Interpreter', 'none');
        box off;
        axis tight
    end


    % plot corresponding ephys stim pulse
    % TODO: determine timedata here, instead of relying on the above
    if currentTrace{1}.header.ephys.ephys.startButton && currentTrace{1}.header.ephys.ephys.stimOnArray(1) == 1;
        try % TODO: get rid of the try 
            traceLength = currentTrace{1}.header.ephys.ephys.traceLength;
            if ~pulseJackerFlag
                params = currentTrace{1}.header.ephys.ephys.pulseParameters{:};
                if isempty(params)
                    try
                        params = currentTrace{1}.header.ephys.ephys.pulseParameters{2}; % probably wrong for 2nd chan
                    catch
                        disp('qvTrace: field pulseParameters is empty.')
                    end
                end
                % here's where the problem is. if there's only one amp channel,
                % there's only one index. but if there are two channels, the
                % first channel's pulse has index 2; index 1 is empty.
            else
                cp = currentTrace{1}.header.pulseJacker.pulseJacker.currentPosition;
    %             np = size(currentTrace.header.pulseJacker.pulseJacker.positions{1}, 2);
    %             params = currentTrace.header.pulseJacker.pulseJacker.pulseDataMap{(cp-1)*np+1}(:);
                channelIndex = 1;
                params = currentTrace{1}.header.pulseJacker.pulseJacker.pulseDataMap{channelIndex, cp + 1};
            end
            pulseTrace = pulseTraceFromParameters(params, traceLength);
            timedata = ((0 : currentTrace{1}.header.ephys.ephys.sampleRate * currentTrace{1}.header.ephys.ephys.traceLength - 1) / ...
                currentTrace{1}.header.ephys.ephys.sampleRate)';
    %         if loopFlag && loopDaqFlag
    % %             pulseTrace = repmat(pulseTrace, loopIterations, 1);
    %             pulseTrace = [];
    %             for n = 1 : loopIterations
    %                 np = size(currentTrace{1}.header.pulseJacker.pulseJacker.positions{1}, 2);
    %                 params = currentTrace{1}.header.pulseJacker.pulseJacker.pulseDataMap{(n-1)*np+1}(:);
    %                 pulseTracelet = pulseTraceFromParameters(params, traceLength);
    %                 pulseTrace = [pulseTrace; pulseTracelet];
    %             end
    %             size(pulseTracelet)
    %             size(pulseTrace)
    %             timedata = ((0 : currentTrace{1}.header.ephys.ephys.sampleRate * currentTrace{1}.header.ephys.ephys.traceLength ...
    %                 * loopIterations - 1) / currentTrace{1}.header.ephys.ephys.sampleRate)';
    %         end
            plotnum = plotnum + 1;
            subplot(Nsub, 1, plotnum);
            hplot = plot(timedata, pulseTrace);
            set(hplot, 'LineWidth', 2, 'Color', 'r');
            xlabel('Time (seconds)');
            ylabel('Amplitude');
            htitle = title(['Pulse: ' params.name]);
            set(htitle, 'Interpreter', 'none');
            box off;
            axis tight
        catch
            disp('pairedAverage: problem making ephys stim pulse.')
        end
    end
end
% ============ EPHYS TRACE, CHANNEL 2 ==================
% TODO: plot multiple channels using FOR loop
% if currentTrace{1}.header.ephys.ephys.startButton && currentTrace{1}.header.ephys.ephys.acqOnArray(1) == 1;
%     
%     ephysdata = zeros(fileLength,length(currentTrace{1}.data.ephys.trace_2));
%     for i = 1:fileLength
%         ephysdata(i,:) = currentTrace{i}.data.ephys.trace_2;
%     end
%     
%     timedata = 0 : size(ephysdata, 2) - 1;
%     timedata = timedata / currentTrace{1}.header.ephys.ephys.sampleRate;
%     plotnum = plotnum + 1;
%     subplot(Nsub, 1, [plotnum plotnum+1]);
%     plotnum = plotnum + 1; % double-high subplot for ephys traces
%     
%     mode = input('Take average (1) or draw each trace (2)? Say 1 or 2: ');
%     switch mode
%         case 1
%             plot(timedata, mean(ephysdata,1));
%         case 2
%             plot(timedata,ephysdata)
%     end
%     xlabel('Time (seconds)');
%     ylabel('Amplitude');
%     traceType = '';
% %     if loopFlag
% %         if loopCpuFlag
% %             traceType = [' -- cpu loop, iteration ' ...
% %                 num2str(loopIterationCounter) ' of ' num2str(loopIterations)];
% %             currentTrace.userFcn.fIAnalysis.familyType = 'cpu loop';
% %         elseif loopDaqFlag
% %             traceType = [' -- daq loop, all ' num2str(loopIterations) ' iterations'];
% %             currentTrace.userFcn.fIAnalysis.familyType = 'daq loop';
% %         end
% %     else
% %         currentTrace.userFcn.fIAnalysis.familyType = 'not a loop';
% %     end
% %     htitle = title([pathstr '\' tracename traceType]);
% %     set(htitle, 'Interpreter', 'none');
%     box off;
%     axis tight
% end
% 
% 
% % plot corresponding ephys stim pulse
% % TODO: determine timedata here, instead of relying on the above
% if currentTrace{1}.header.ephys.ephys.startButton && currentTrace{1}.header.ephys.ephys.stimOnArray(1) == 1;
%     try % TODO: get rid of the try 
%         traceLength = currentTrace{1}.header.ephys.ephys.traceLength;
%         if ~pulseJackerFlag
%             params = currentTrace{1}.header.ephys.ephys.pulseParameters{:};
%             if isempty(params)
%                 try
%                     params = currentTrace{1}.header.ephys.ephys.pulseParameters{1}; % probably wrong for 2nd chan
%                 catch
%                     disp('qvTrace: field pulseParameters is empty.')
%                 end
%             end
%             % here's where the problem is. if there's only one amp channel,
%             % there's only one index. but if there are two channels, the
%             % first channel's pulse has index 2; index 1 is empty.
%         else
%             cp = currentTrace{1}.header.pulseJacker.pulseJacker.currentPosition;
% %             np = size(currentTrace.header.pulseJacker.pulseJacker.positions{1}, 2);
% %             params = currentTrace.header.pulseJacker.pulseJacker.pulseDataMap{(cp-1)*np+1}(:);
%             channelIndex = 2;
%             params = currentTrace{1}.header.pulseJacker.pulseJacker.pulseDataMap{channelIndex, cp + 1};
%         end
%         pulseTrace = pulseTraceFromParameters(params, traceLength);
%         timedata = ((0 : currentTrace{1}.header.ephys.ephys.sampleRate * currentTrace{1}.header.ephys.ephys.traceLength - 1) / ...
%             currentTrace{1}.header.ephys.ephys.sampleRate)';
% %         if loopFlag && loopDaqFlag
% % %             pulseTrace = repmat(pulseTrace, loopIterations, 1);
% %             pulseTrace = [];
% %             for n = 1 : loopIterations
% %                 np = size(currentTrace{1}.header.pulseJacker.pulseJacker.positions{1}, 2);
% %                 params = currentTrace{1}.header.pulseJacker.pulseJacker.pulseDataMap{(n-1)*np+1}(:);
% %                 pulseTracelet = pulseTraceFromParameters(params, traceLength);
% %                 pulseTrace = [pulseTrace; pulseTracelet];
% %             end
% %             size(pulseTracelet)
% %             size(pulseTrace)
% %             timedata = ((0 : currentTrace{1}.header.ephys.ephys.sampleRate * currentTrace{1}.header.ephys.ephys.traceLength ...
% %                 * loopIterations - 1) / currentTrace{1}.header.ephys.ephys.sampleRate)';
% %         end
%         plotnum = plotnum + 1;
%         subplot(Nsub, 1, plotnum);
%         hplot = plot(timedata, pulseTrace);
%         set(hplot, 'LineWidth', 2, 'Color', 'r');
%         xlabel('Time (seconds)');
%         ylabel('Amplitude');
%         htitle = title(['Pulse: ' params.name]);
%         set(htitle, 'Interpreter', 'none');
%         box off;
%         axis tight
%     catch
%         disp('pairedAverage: problem making ephys stim pulse.')
%     end
% end
title(fileLength)